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I. INTRODUCTION 



Supersymmetry (SUSY) is playing an increasingly relevant and unifying role in high 
energy physics. Prom a purely theoretical point of view, SUSY is required for consistency 
and finiteness in superstring theory; compactification and SUSY breaking mechanisms are 
then needed in order to produce a low energy four-dimensional effective action with a residual 
A'" = 1 SUSY. This constraint comes from the phenomenological side where the most popular 
current extensions of the Standard Model are actually based on SUSY for at least two 
reasons. First, supersymmetric Grand Unification theories are quite successful in predicting 
SU{3) X SU{2) X U{1) gauge couplings unification [1], a fact that can be considered as the 
main phenomenological motivation for SUSY [2]. Moreover, supersymmetric models solve 
in a natural way the hierarchy problem [3] of matching the electroweak and GUT scales 
without being spoiled by huge radiative corrections to Higgs masses. 

However, many features of this scenario still need some clarification. Indeed, N = 1 
SUSY is expected to be exact at the GUT scale around 10^^ GeV, but must be broken in 
the TeV region in order to account for the asymmetric mass textures that are currently 
observed. In particular, this will be true if some superpartner with a mass of a few TeV will 
be observed in the future LHC or Linear Collider experiments. In the model independent 
analysis, the source of breaking is usually parametrized by a complete set of soft terms whose 
origin remains however rather unexplained. 

In several approaches, it is due to some kind of spontaneous breaking of SUSY in a 
hidden sector and communicated to the MSSM particles producing the soft terms. As with 
every dynamical symmetry breaking, non-perturbative techniques must be exploited and 
the lattice regularization and renormalization programme is of course one of the main lines 
of research. Indeed, the simultaneous introduction of infrared and ultraviolet cutoffs allows 
for calculations, like strong-coupling expansions, that are quite complementary to the usual 
weak-coupling perturbative analysis. 

Beside this, when any known analytical treatment fails, lattice models can also be studied 
by direct simulations that can provide, in principle, accurate numerical measurements. 

In this paper, we address the problem of spontaneous supersymmetry breaking (S^B) 
in a simple, but interesting, theoretical laboratory that is the class of Wess-Zumino (WZ) 
two-dimensional models of chiral superfields with no vector multiplets. Preliminary results 
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on this subject can be found in [4]. Related studies of the two dimensional Wess-Zumino 
model can be found in [5]. 

Despite their simplicity, these systems are nontrivial because in two dimensions super- 
symmetry is not strong enough to predict the exact pattern of breaking, a situation that 
must be compared with the four dimensional case where WZ models are expected to break 
supersymmetry if and only if they do at tree level. 

Unfortunately, as we shall discuss, lattice strong-coupling expansions provide useful in- 
sights, but are unable to reliably predict the physics of the continuum theory and one must 
resort to a numerical analysis. 

Since S^B is closely related to the symmetry properties of the ground state, it appears to 
be quite reasonable to adopt some kind of Hamiltonian formulation. Moreover, we will see 
in the following that, if we wish to preserve a SUSY subalgebra, a conserved Hamiltonian 
is crucial. However, the traditional algorithms for simulation of lattice field theories are 
based on the Lagrangian formulation [6]. The main reason is the immediate probabilistic 
interpretation of the partition function, at least for bosonic systems not suffering from a 
sign or phase problem; this leads to a host of Monte Carlo algorithms, some of which are 
extremely efficient. Of course, alternatives based on a more direct Hamiltonian formalism 
do exist [7], but they are certainly less exploited in high energy physics where emphasis is 
on Lagrangian symmetries, in particular Poincare invariance. 

On the other hand, Hamiltonian methods has been used in Supersymmetric Discretized 
Light-Cone Quantization (SDLCQ) [8] and also are widely exploited in non relativistic con- 
texts [9] where the properties of the ground state are typically the simplest and first object 
of investigation. Moreover, these techniques interlace the brute force numerical calculations 
with analytical or physical insights about the structure of the ground state wave function, 
a feature that is quite welcome in the study of S^B where we expect major changes to show 
up at the phase transition. 

Another important feature of our study concerns the fact that fermions, needed in su- 
persymmetric models, are the major source of complications in the current approaches to 
lattice simulations. In the Lagrangian approach quantum expectation values are computed 
by summing over histories of the classical fields, following Feynman's ideas. In the case of 
fermions, these are Grassmann valued classical fields that cannot be analyzed by direct nu- 
merical methods. The typical solution amounts to integrate them out and study the resulting 
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non-local bosonic model [10]. This can be nontrivial for a generic model and a recent de- 
tailed account of this problem and whether it can be formulate successfully supersymmetric 
theories on the lattice can be found in [11]. 

Instead, in the Hamiltonian approach, what is relevant is the algebra of the fields and 
their conjugate momenta. From this point of view, fermions and bosons differ just by the 
replacement of commutators by anticommutators and also by the structure of the state space, 
finite dimensional for fermions in finite volume, infinite dimensional in the bosonic case. 
Apparently, in the Hamiltonian approach, there is much more symmetry in the treatment 
of fermions and bosons than in the Lagrangian approach. 

Looking at the details of the simulation techniques, however, problems arise with Hamil- 
tonian fermions due to the well known hard sign problem [12]. Roughly speaking, fermion 
exchanges introduce problematic and unavoidable signs that often break in a substantial way 
the probabilistic interpretation of quantum expectation values required to build a numerical 
stochastic algorithm. This deep problem is somewhat milder in 1 -|- 1 dimensions where 
specific equivalences between fermionic and bosonic fields can be estabhshed [13]. Also, the 
topology of fermion dynamics is the simplest possible and helps in taming the sign prob- 
lem. Actually, for several fermion models in 1 + 1 dimensions arising in Solid State theory, 
like, e.g., Hubbard-like models, algorithms can be devised with no sign problem and good 
efficiency as well as scahng properties [15]. 

The detailed plan of the paper is the following: in Sec. (II) we present the model and its 
lattice Hamiltonian; in Sec. (HI) we compute the first nontrivial order of the strong-coupling 
expansion of the ground state energy; in Sec. (IV) we discuss the Renormahzation Group 
trajectories along which a continuum limit can be reached. In Sec. (V) we describe a Green 
Function Monte Carlo algorithm. Finally, Sec. (VI) is devoted to present our numerical 
results. 

II. THE iV = 1 WESS-ZUMINO MODEL IN 1 -|- 1 DIMENSIONS 
A. Definition of the model and patterns of SUSY breaking 

Let us consider the most general SUSY algebra in two dimensions. The generators are 
split into fermionic and bosonic ones. The algebra with N left-handed fermionic generators 
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{Ql}a=i,...,n and N right-handed fermionic generators {Qr} a=i,...,n is denoted by {N,N). 
The bosonic generators are the components of the two-momentum (P°, P^) and a set of 
central charges T"^^. The non trivial part of the algebra is 

{QlQ^n} = 6^^{P' + P'), 

In the left-right symmetric case with {N, N) = (1, 1), we denote 

Q^,2^Q]^±Ql, (2.1) 

and find 

{Qa, Qb} = 2{H1 + Pa^ + Ta^),,, (2.2) 

where cr* are the Pauli matrices, (P°,P^) = {H,P) and T = T^^. The minimal realization 
of this algebra requires a single real chiral multiplet with a real scalar component (p and a 
Majorana fermion with components ■0^'^. The explicit form of the supercharges is 



Qm = / 



dx 



(2.3) 



where p{x) is the momentum operator conjugate to (p{x). The central charge corresponds 
to a topological quantum number [16] 

r = /d.gn*'). (2.4) 

As usual with supersymmetric models, the structure of the Hamiltonian H guarantees that 
the energy eigenstates have £^ > because 

H=\{Qi + Ql). (2.5) 

Invariant states annihilated by both Qi coincide with the zero energy states and arc thus 
supersymmetric ground states; they must lie in the topologically trivial sector. 

The problem of predicting the pattern of S^B is not easy. In principle, the form of 
V{(f) is enough to determine whether supersymmetry is broken or not. At least at tree 
level, one easily proves that supersymmetry is broken if and only if V has no zeros. In 
two dimensions, however, these conclusion is generally false due to radiative corrections. 



An analytic non-perturbative tool that can help in the analysis is the Witten index defined 
as [23] 

/ = Tr(-l)^, (2.6) 

where F is the fermion number. Since super symmetry is not explicitly broken, contributions 
from positive-energy bosonic and fermionic states cancel and 

/ = - nLo- (2.7) 

In finite volume, / is invariant under changes in V{ip) that do not modify its asymptotic 
behaviour. In particular, it can be computed at weak coupling where each zero of V{i~p) is 
associated to a perturbative zero energy state. Thus, if V{i~p) has an odd number of zeroes, 
we find / 7^ and supersymmetry is unbroken. If, on the other hand, V{i~p) has an even 
number of zeroes, the associated perturbative vacua can contribute / with opposite signs 
and, when / = 0, we cannot say anything. In particular, a nontrivial set of perturbative zero 
energy states with / = can receive instanton corrections due to tunneling lifting them to 
positive energies breaking supersymmetry while leaving 7 = 0. In such cases, the behaviour 
of the tunneling rate with increasing volume is crucial in answering the question of breaking. 

An interesting example of this complicated scenario is discussed in Appendix A of Ref. 
[23]. We quickly review the analysis since it will be important in the interpretation of our 
results. When V{(p) = A((/9^ + a^), the action of the WZ model is 

S^Jd'x (^\{dcpy + Iv^T • di: - ^X'iip' + - ^A^V^V') ■ (2.8) 

For large positive the index is zero because there are no zero-energy states. Due to a 
special conjugation symmetry valid for this model in finite volume, the pattern of breaking is 
invariant under —a^. This means that for negative a^, the two zeroes of V are bosonic 
and fermionic and (finite volume) tunneling lifts their energy to a positive value. However, 
in infinite volume and large negative a^, the narrow minimum of the potential is protected 
from radiative corrections and generates an expectation value {(p) ^ signaling the SSB of 
the Z2 symmmetry >^ —99, ■?/' 75?/'. The fermion becomes massive and supersymmetry 
must be unbroken due to the absence of a massless Goldstino. 

The above discussion illustrates that an alternative non-perturbative analysis of the mod- 
els with 7 = is certainly welcome. In the following, we shall put the model on a space-time 
lattice in order to exploit explicit numerical simulations as well as analytical strong-coupling 
expansions. 
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B. Lattice Version of the Model 



On the lattice it is impossible to maintain the full SUSY algebra and it is important to 
understand what can be said by looking at subalgebras. If we consider one supercharge only, 
for instance Qi, and find a state with Qi\s) — 0, we cannot say that it is a zero energy state 
unless we know that it is in the T = sector. On the other hand, if no states with Qi\s) = 
are found in any topological sector, then supersymmetry is certainly broken. 

Thus, even if we forget Q2, we can choose as a clear-cut signal of supersymmetry dynam- 
ically breaking the lowest eigenvalue of the operator Qf: if it is positive, we have breaking. 

The SUSY algebra (2.2) involves exphcitely the generators of space and time infinitesimal 
translations and is spoiled on the lattice. In the Lagrangian approach, both space and time 
are discrete and SUSY is completely broken. In the Hamiltonian formulation, time remains 
continuous and the D = 2 algebra is reduced to = 1 and not totally lost. The full 
two-dimensional algebra as well as Lorentz invariance are expected to be recovered in the 
continuum limit. 

A lattice version of the above model has been previously studied in Refs. [17, 18]. A 
similar approach to Wess-Zumino models with N — 2 supersymmetry is discussed in Ref. 
[19], and numerical investigations are reported in Ref. [20]. On each site of a spatial lattice 
with L sites, we define a real scalar field together with its conjugate momentum such 
that [pn, ^pm] = —i^n^m- The associated fermion is a Majorana fermion with a = 1, 2 
and {'0a,n,'06,m} = 5a,b5n,m , V'ln = V'a.n- The fcrmiouic charge 

Q = E 

n=l 

with arbitrary real function V{(p), (called prepotential in the following) can be used to define 
a semi-positive definite lattice Hamiltonian 

H - (2.9) 

This Hamiltonian includes the central charge contribution in the form of a term 

Eny^n) ^"^^:^"-\ (2.10) 
n ^ 

that is precisely a discretization of T. Eigenstates of B. are divided into invariant Q-singlets 
with zero energy and Q-doublets with degenerate positive energy. H describes an interacting 
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model, symmetric with respect to Q and this symmetry is respected by the spectrum if and 
only if the ground state energy is positive. We stress again that Q symmetry breaking 
implies breaking of the full 2 dimensional supersymmetry, whereas Q symmetry does not 
imply (in a generic topological sector) 2D SUSY. 

We remind that, on the lattice, spontaneous supersymmetry breaking can occur even for 
finite lattice size L, because tunneling among degenerate vacua connected by Q is forbidden 
by fermion number conservation. 

To write in a more famihar form, we follow Ref. [18] and replace the two Majorana 
fermion operators with a single Dirac operator x satisfying canonical anticommutation rules, 

i-^-) {Xn; Xm} 0) {Xnj Xml ^n,m- 

A,n = ^"^^."' (Xl + iXn), V'2,n = ^^^^^ (xl ' ^Xn)- (2-11) 

The Hamiltonian takes then the form 

-\{xiXn+i + h.c.) + {-lTV'{<p^) {xiXn - ^)} (2.12) 

and describes canonical bosonic and fermionic fields with standard kinetic energies and a 
Yukawa coupling. 

This Hamiltonian conserves the total fermion number 

Nf^Y^xiXn, (2.13) 

n 

and can be examined in each sector with definite Nf separately. For reasons that will be 
understood later, we shall also consider open boundary conditions and restrict the lattice 
size L to be even. These are constraints that do not affect the physics of the model in the 
continuum, but will be very welcome by the algorithm we are going to apply. 

The simplest way to analyze the pattern of supersymmetry breaking for a given V{ip) is 
to compute the ground state energy Eq. As wc mentioned, on the lattice, we can perform 
such a computation in a non-perturbative way by strong coupling or numerical simulations. 
However, before discussing these items, we want to stress some identities that can be used 
together with Eq to get information on the symmetry of the ground state. 
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C. Lattice Ward identities 



If the vacuum |0) is supersymmetric, Q\0) — and for each operator X we have 

(0|{g,X}|0) = 0. (2.14) 

In particular, taking 

X^J2FMi^2,n, (2.15) 

n 

we obtain 



{o\j:\F{<pn) 



'^"^^ ^ '^""^ + ^(^n)] + F'{<p^){-ir{xixn - 1/2)} |0) = 0. (2.16) 

A basis of Ward Identities is thus obtained by considering F{ip) — ip'^. For instance, on 
an even lattice with open boundary conditions we find for n = 1 the relation 

(0| Y: {^nVi^n) + i-lTxiXn} |0) = 0. (2.17) 
n 

The case F{(p) — constant is also interesting. It leads to 

{0\J2VM\0) = 0. (2.18) 



III. STRONG COUPLING ANALYSIS OF SUSY BREAKING 



Let us start from the supersymmetry charge 

Following Ref . [19] , we define the strong-couphng hmit by 

y((^) — Ai-(°)(A(^), 

A— >oo 



perform the canonical transformation 



(0) 



A 



and rescale the energies by A^; dropping the index from and p, the result is 

(mi - "Pi-^Wx 



2A2 



pf + V\ipi) + 2iV'{^i)4^l;f 



A2 



4A4 



Introducing the Dirac fields Xh xl [18]; cf. Eq. (2.11), we obtain 

H^'^ - E [IpI + Iv'i^i) + (-l)VV/)(xb - 1/2) 



^^'^-^E(^Hi-^z-i)^ 

where we denote by n; = 0, 1 the eigenvalue of xlxi- 



A. Leading order 

To leading order in 1/A, the Hamiltonian is factorized in a supersymmctric quantum me- 
chanics for each site; adopting an explicit representation, we can write the one-site Hamil- 
tonian as 



(3.1) 



(in the occupation number representation the basis chosen now is (n=0, n—1) for odd sites 
and (n=l,n=0) for even sites). This Hamiltonian has a, N — 2 supersymmetry [21]: 



{Qh Qj} = 5ijH, Qi = \ + (y2V{x)] , Q2 = \ [cr2P - criV{x) 
The conditions for a supersymmctric ground state Qiipo — reduce to 

i/j'oix) = a3V{x)'ilJo{x). 



(3.2) 



(3.3) 



For polynomial V, supersymmetry is unbroken if and only if it is possible to find a normal- 
izable solution to Eq. (3.3), which happens if the degree of F is odd [21]. 
We can write the time-independent Schrodinger equation as 

V'" + [2E - V^{x) T V'{x)) V' = 0; 

denoting the eigenfunctions for the two signs by and their energies by we have 

+ (24 - V\x) T V'{x)) = 0. (3.4) 

Supersymmetry implies that, for E 0, states are paired in boson-fermion doublets. 
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FIG. 1: The ground-state energy eq = and the overlap tjq as functions of Aq for the quadratic 
prepotential V = ip'^ + Xq. 

We remark that this conclusion is in strong disagreement with the continuum (or weak- 
coupling lattice) analysis where the relevant feature of V is the existence of zeros. 

For g > 1, ^m(^) cannot be computed exactly (excluding the cases e"^ = 0); it 

is however easy to compute then numerically to high accuracy, using, e.g., the Numerov 
method [22]. An example is shown in Fig. 1. 

In the following analysis, we shall have to tell between V with odd or even leading power 
of (f. 

1. Odd q 

For odd q, we either have Eq = or Eq = 0; let us assume Eq = 0: ipo is the supersym- 
metric ground state satisfying Eq. (3.3); all the other states appear in pairs: e^+i = e^. 
Notice that the ground state is bosonic for even sites, fermionic for odd sites (the opposite 
holds if 4 = 0). 
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A strong argument against supersymmetry breaking is given by the Witten index 
Iw = Tr{—l)^f [23]; in the strong-couphng hmit, we clearly have Iw — ±1; since Iw 
implies unbroken SUSY, and Iw is invariant under regular perturbations (cf. Sect. Ill A 3), 
supersymmetry can never be broken for odd q, not even in the L ^ oo limit. 

A simple check of this result can be given explicitely when V is linear and is discussed in 
details in App. (A). 



2. Even q 

For even g, we have £^ — ^m'^ m = 0, this corresponds to a degenerate ground state 
with broken supersymmetry [sq = > 0). The phases of the normalized states \'ipn) 
be chosen in order to satisfy 



/2^|V^+) = [zp + n(^)]|V^-); 



(3.5) 
(3.6) 



Introducing the notation 



(^o'lV'o") = ^0, 



we can prove the important relations 



fracl\/2s^r]o. 



(3.7) 
(3.8) 



?7o can be computed numerically from ip^(Lp), cf. Fig. 1. The proof of Eq. (3.7) is very simple: 
just take the scalar product of Eq. (3.5) with {iPqI and of Eq. (3.5) with {iPq\, and observe 
that (p)± — 0. The proof of Eq. (3.8) is also immediate: 

Several simplifications can be exploited when V{ip) is even. For an asymptotically positive 
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polynomial V{(p) with degree g > 2 it is easy to show that ^ 

where / is the hermitian parity inversion operator 

{<p\m = {-<p\^|;) 
Then, the eigenstates can be characterized by the single equation 

where 

It is easy in this case to obtain generalized relations like the previous ones. Let us consider 
the equation 

V2i:{lPm\mn) = {-iniPm\Iiip + V)I\lPn) = (-1)" + 1^) IV'n) 

Taking the hermitian conjugate and exchanging m and n we find the two equations 

V^(V'n|/|^™) = + (3.9) 

V2^{ll^nmm) = {-ir{llJn\Hp + V)\xPm) (3-10) 



therefore 



or (exploiting parity) 

mv{^)\^^) = ^iv^+v^i-ir^'^mii^^) (3.11) 

In a similar way we can compute 

V^(V'mklV'n) = {-lT{^|JmH^p + V)I\'^|Jn) 

Taking the hermitian conjugate and exchanging m and n we find the two equations 

A/2i;(V'n|<^|V'm) = (-l)"(^n|(-^P - V)'^I\^^) (3.12) 

V2iZ{^n\^\^m) = {-in^l^nMip + V^)/|^m) (3.13) 



^ In fact, from their definition, one can see that ^^{(f) have the same sign when (p +oo. Since they are 
related by a parity transformation, their relative phase is fixed by the number of nodes. 
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summing the two equations 
and (exploiting parity) 

Of course, for n = m = 0, Eqs. (3.11), (3.14) agree with the previous general results. 
3. On the convergence of the perturbative expansion 

The Rayleigh-Schrodinger perturbation theory of a Hamiltonian of the form H — Hq + 
PHi is regular (i.e., it has a finite radius of convergence in /3) if [24] 

< a||i/o*|| + (3.15) 

uniformly for all state vectors in our case, Eq. (3.15) clearly holds, for both H^"^^ and 
H^'^\ except for the trivial cases q < 1. 

B. First order 

At first order (subleading) in the strong-coupling expansion we consider again the two 
cases of even or odd q. 

1. Odd q 

In the case of unbroken supersymmetry (odd q) , the subleading correction to the ground- 
state energy in the strong-coupling expansion is zero: the fermionic contribution iil)l_^-^il)f + 
ii^i+ii^i has clearly zero diagonal matrix elements, and the bosonic contribution (pi^iV{(pi) — 
(piV{ipi+i) is zero because it factorizes into {<p){V) — {(p){V); strictly speaking, this is true 
for periodic and free boundary conditions, but it could be false, e.g., for fixed boundary 
conditions. 
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2. Even q 



Due to the structure of the Hamiltonian, it is convenient to describe states in the mixed 
form 

Yl V'ni,...,ni(<^l,---,</'L)|ni,...,nL) (3.16) 
ni,...,nL 

where ipni,...,nL{fiT--ifL) is a wave function depending on the bosonic degrees of free- 
dom and \ni, . . . jUl) is the fermionic component of the state defined in terms of the state 
annihilated by all x 

X^\0) = 0, (3.17) 
according to the canonical ordering of the Fermi operators: 

\ni,...,nL)^{x\r---ixir''\Q)- (3.18) 

Of course = 0, 1 and |ni, . . . , ul) describes a state with rii fermions at site i. In the case 
of broken supersymmetry (even q) , the subspace B of lowest leading-order energy is spanned 
by the states 

|n) = V'o'((/Pi) • • •<M</'l) Im, . . . (3.19) 

where ai = (—1)'+'*' and ip^^ = ip^. We have adopted open boundary conditions, corre- 
sponding in our notations to setting (po — (fiL+i = and ipQ = ipl+i = (and therefore 
Xo = Xl+1 = 0). 

Since the number of fermions "^^i is conserved, we can impose an additional constraint 
on B and define 

BN = {\n),Y.'ni = N}, B = Bo®...®Bl. 
I 

We will now prove that, for even L, the ground state is doubly degenerate and lies in the 
sectors with N = L/2, L/2 - 1. 

At first order, we have to diagonalize the operator H^'^^ over Bn- Let us spht 

= //g) + (3.20) 

H^B = ^i:^(¥';)(mi-<^/-i) (3.21) 
^ 1=1 

hP = -I Eixlxi+i + xUiXi) (3.22) 
^ 1=1 



15 



Since 



we have 



{n'\H^i\n) ^ --A,r.' [(-ir + (-i; 



where we have exploited 



V2£o 



that holds for even V. Since n = 0, 1 we can use (—1)" — l — 2n and write 



The matrix elements of H^I^ are 



(3.23) 
(3.24) 
(3.25) 

(3.26) 
(3.27) 



where /iri,Ti' = 1 if n and n' are connected by Hp^ (i.e. a hopping of one fermion) and 
otherwise. 

Thus, we can hide the bosonic wave functions and write an effective perturbation acting 
on purely fermionic states as 



(3.28) 



where 1 is the identity operator and 



-1 \t-3\ = l 
1 i — j — 1 or L 
otherwise 



(3.29) 



Since H^^ is quadratic, it is convenient to change operator basis. Let v^^^ be the p-th 



eigenvector of M with eigenvalue A*^^^: 



12 -S, 



p,L 



A(-) = -2r-^^ 



L 

:cos- 



sm 



pn 
~L 



(3.30) 
(3.31) 



They define a (real) unitary matrix 



p=\ i=l 
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We can replace the operators Xi by the operators Oj defined by 

L L 



p=l i=l 



with 



The new form of H^q is 



pq 



HS-i\j:>^^H<^p-^] (3-32) 

ip=i 



The eigenvalues X^^\ . . . , X^^/'^~^^ are negative and X^^/^^ — 0; there are thus two degenerate 
ground states with L/2 — 1 and L/2 fermions. This is required by supersymmetry: since 
H^"^^ restricted to B commutes with Q^°^ , all the states must be paired in doublets with N 
differing by 1. The ground state with L/2 fermions is 

= Vo'M ■ ■■roH<PL) a{ ■ • -a^/^IO > (3.33) 
To conclude, the shift of the ground state energy due to the perturbation is 

^..|L_,f:.o.M._fot^ (3.34, 



n=l 



In the L — > oo limit we have 



^^-!i + 0{l/L) (3.35) 



In summary, the first order perturbation in the strong-coupling expansion removes the large 
degeneracy of the ground state and determines a doublet of eigenstates with L/2 — 1 and 
L/2 fermions with minimum energy 

^ = s.-l'4^0(4r,l) (3.36) 



L ^ A2 TT KX^L' X\ 

A similar calculation at first order for {(pk) and {(pk^i)c is discussed in App. (B). The second- 
order correction to the ground state energy can also be computed with a reasonable effort 
and the result is described in App. (C). However, we remark that the results drawn from 
the first order corrections are not qualitatively changed. 
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C. Discussion 



Prom the analysis of the strong-couphng expansion of the model we can draw the following 
conclusion. For polynomial V{(fi), the relevant parameter is just its degree q. 

For odd q, the strong coupling analysis and the tree-level results agree and supersymmetry 
is expected to be unbroken. This conclusion gains further support from the nonvanishing 
value of the Witten index at strong coupling. 

For even q, in strong coupling, the ground state (at least in the sector with half filling) 
has a positive energy density also for L — > oo and supersymmetry appears to be broken. 
Of course, this can be in disagreement with weak coupling. A specific case that we shall 
analyze numerically in great detail is 

V{(f) = X2<f^ + Xq. (3.37) 

For Ao < 0, weak coupling predicts unbroken SUSY, whereas the strong coupling prediction 
gives broken SUSY for all Aq. For this specific model, as we already discussed, the strong 
coupling analysis agrees with Ref. [23] in the sense that it reproduces the continuum physics 
in finite volume. 

For large expansion parameter, the strong coupling results can be compared with explicit 
simulations (that we shall fully discuss in Sec. V). Let us consider for instance the quadratic 
model with Aq = on a lattice with L — 22 spatial sites. In Fig. 2, we show the expectation 
value {(fn) computed at A2 = 2. The agreement is quite good apart from the points on 
the border where the convergence seems to be slower. To check the vahdity of the strong 
coupling expansion at smaller couplings, we show in Fig. 3 the ground state energy from 
MC simulation compared with the first and second order strong coupling expansion. The 
scaled variables on the plot axes are discussed in App. (C). The second order gives better 
results at large values of the expansion parameter, but is unreliable below A2 — 0.35. 

In the next Section, we shall see that the continuum limit is in the region of small A2. 
Thus, for even q, it seems difficult to gain additional insight from strong coupling and some 
kind of transition can happen as the continuum limit is reached. For this reason, a full 
simulation of the model appears to be the only way to answer the question of symmetry 
breaking. 
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FIG. 2: Comparison between strong coupling and MC simulation for the expectation value 
in the quadratic model with V{ip) = 2ip^ on a lattice with L = 22 spatial sites. 

IV. RENORMALIZATION GROUP TRAJECTORIES 



The action of the WZ model is 

•1 



At the perturbative level, this is a superrenormalizable field theory that can be made finite 
by a renormalization of V{ip). In the minimal subtraction scheme the renormalized potential 
is obtained by solving the heat equation [25] 



1 



(4.1) 



where /i is the dimensional regularization scale. A dependence on /i is thus introduced in 
the coefficients of the various monomials appearing in the tree level V{(f). For the specific 
case of V{(p) = A2V?^ + Aq, we find that A2 is scale-independent and 

Ao(/i) = Ao(Aio) - ^ log —. 

2vr /io 
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FIG. 3: Comparison between strong coupling and MC simulation for the ground state energy in 
the quadratic model with V{ip) = A2 on a lattice with L = 22 spatial sites. 



On the lattice, let us denote by a hat the adimensional lattice coupling constants and by 
the label "ph" the physical ones, fixed and with dimension 1. The above result leads to 



A. 



aAg'^^ Ao--^ log aM. 

ZTT 



(4.2) 
(4.3) 



The way we read these equations is: at one loop and for small enough a, the physical Aq is 
obtained by compensating Aq by the effect of the one-loop diagrams. These are computed 
with the UV cutoff a and with the IR cutoff given by the (dimension 1) mass M of the 
virtual particles in the loop. 

The first equation allows to replace a by A2 everywhere and we get 



A2 
Ao 



Xf A2, 
A.3i + ^log A2-, 



(4.4) 
(4.5) 
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This seems to show that the continuum hmit can be reached with A2 — > and 

i^'^~°A+^logA2 (4.6) 
where A contains the ratio A2^/Ao'^ and the details of the physical mass generation. 



V. SIMULATION ALGORITHM 

A. Green Function Monte Carlo: general considerations 

In this Section, we review the Green Function Monte Carlo approach to the study of the 
ground state of a general quantum model. To this aim, we consider the simple case of + 1 
dimensional quantum mechanics in order to illustrate the basic ideas without unnecessary 
details hiding the essential features of the algorithm. 

For a canonical spinless quantum particle, the Hamiltonian is 

H^^p^ + V{q), [qi,pj] = tSij. (5.1) 

The ground state |^o) of H can be projected out of any initial state |i) with non zero overlap 
{^oli) 7^ 0. The projection is performed by applying the evolution semigroup {exp{—tH)}t>o 
and going to asymptotically large times. 

Focusing on the ground state energy £'0, this procedure leads to the following simple 
formula 

the final state |/) is in principle arbitrary, as long as it is not orthogonal to |^o; in practice, 
it must be chosen with care, to avoid numerical instability. 

The vacuum expectation value of a generic observable O can be computed as 

<*o|Q|*»>=.,!'j? J(;|e-."o«|.)' - 

this procedure is known as forward walking. 

To translate the above formula into a stable numerical algorithm, it is necessary to find 
a basis such that the Hamiltonian H has non positive off-diagonal matrix elements. By the 
way, this is true for the Hamiltonian (5.1) in the basis {\q)} of position eigenstates. If such 
a basis is found, it is possible to identify matrix elements of e""*^ as probability transitions 
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defining a Markov random process in the state space. For instance, in the simple case when 
I/) is chosen to be a zero momentum state, p\f) — 0, we have (Feynman-Kac formula) 

Eo= lim lim J- , (5.4) 

where 'Dqi^r) is the Wiener measure. 

The probabilistic interpretation of the above equation is as follows: £"0 (as well as other 
observables) can be computed by taking the average over weighted walkers which diffuse 
according to the Wiener process. Each path is weighted by the following functional of the 
trajectory 

WlqM] = exp - r V(q(T))dT. (5.5) 
Jo 

In the numerical implementation, an estimate of £'0 is obtained by computing 



E (V(qt)Wt 

lim ^ !1\ ^ (5.6) 



t— >+oo 



where qt is a numerical discretization of the Wiener process, Wf its associated weight com- 
puted by properly approximating Eq. (5.5) and, finally, E (■) denotes the average with respect 
to the realizations of qt. In practice, after the choice of a particular approximation qt, one 
works with a large number K of walkers and extrapolates numerically to i^' — > 00. The 
control of the approximations involved in this strategy requires some discussion that we 
defer to the Section devoted to results. 

A point that is worth mentioning regards the possibility of introducing a guidance in the 
walkers diffusion. To improve the convergence to ground state it is customary to define the 
unitarily equivalent Hamiltonian 

H^e^He-^ ^^p'^ + ip-F + V{q), (5.7) 

where 5" is an arbitrary (real) function and 

F = VS, (5.8) 

It can be shown that the derivation of expressions like (5.4) can be easily generahzed to this 
case and the required modifications can be summarized as: (i) the potential V is replaced 
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by with V, (ii) the Wiener process is replaced by a deformed process guided by the drift F. 
In the following, we shall call Importance Sampling the trick of exploiting a non zero S. 

In the following Sections, we describe in full details the algorithm for the model un- 
der study considering first the bosonic and fermionic sectors separately and finally the full 
Hamiltonian. 

B. Bosonic sector 

The bosonic sector of the lattice model is a canonical quantum mechanical one with many 
degrees of freedom. The algorithm is the one described in the previous Section. Given the 

transformed Hamiltonian H of Eq. (5.7), we write 

eM-eHs) = exp exp e'-^-^exp exp + 0{s'). (5.9) 

The function V depends on the bosonic state, i.e. the set of values of the scalar fields that 
we collectively denote by Q. 

The update rule for the weighted walker {Q, W) is built, step by step, following the 
approximate operator factorization (5.9) and reads (see [26] for similar calculations in the 
solution of the Langevin equation) 

{Q\W')^{Q"',W"'), (5.10) 
where Q'" and W" are built according to 

W"^W'exp(^-'-ViQ')y (5.11) 

Q" ^Q'+ e', 

z^Q'+'-FiQ'), 
Q" = Q' + eF{z), 

Q'" = Q" + e", 

and ^' ^ C," are independent sets of Gaussian random numbers. In the above update, the 
integration of the equations of motion associated with the driving force F has been solved 
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at second order. In the end a systematic error 0{e^) with respect to the evolution time has 
been introduced. 

An estimate of the energy in the bosonic sector is obtained by taking the weighted average 
of V over several realizations of the walker path 



C. Fermionic sector 

In the fermionic sector, the spirit of the algorithm is the same, but there are important 
technical differences that we want to emphasize. At fixed scalar fields configuration, the 
remaining state space is purely fermionic and, on a finite lattice, it is both discrete and 
finite dimensional. 

To simplify notation, in this Section we denote H = H}, . The Hamiltonian can be 
thought as a large sparse matrix H — ||i?ss'|| with s and s' denoting fermionic states. We 
now show that a similar construction like the one exploited with Hb can be repeated here. 
The Gaussian random noise that was the building block in the simulation of the Wiener 
process is replaced here by a discrete jump process. 

Again, the problem is that of giving a probabilistic representation for the evolution semi- 
group Q = {e^*'^}t>o- To this aim, we define a Markov process that describes diffusion in 
the discrete state space and also provide a rule for the update of a walker weight. We finally 
show that suitable averages over weighted walkers reconstruct the evolution governed by H 
and project a given initial state onto the ground state. 

For each pair s, s' in the state space S such that s ^ s' and Hg's 7^ we define Vg's — 
—Hs's- We assume that all Ts's > (no sign problem) and build a 5"- valued Markov stochastic 
process St by identifying Ts's as the rate for the transition s s'. Hence, the average 
occupation Ps{t) — E (^^ , with E(-) denoting the average with respect to St, obeys the 



Related to st, we also define the real valued stochastic process Wt — exp [— JgCOst dtj, 
with cus — ^s'es^s's- It can be shown that the weighted expectation value ipsit) — 
'E{5s,stWt) reconstructs Q,: 




(5.12) 



Master Equation P,(/3) = Es'M^ss'Ps' - ^s'sPs)- 



d 
~dt 



Hss'l/Js'it), 



s' 
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with ■05 (0) = Prob(so = s). Matrix elements of ft can be identified with certain expectation 
values. In particular, the ground state energy Eq (in the purely fermionic sector) can be 
obtained by 

E„ = (5.13) 

The actual construction of the process is straightforward. A realization of St is a piece-wise 
constant map R — > 5" with isolated jumps at times t — to, ti, . . ., with < < ^2 < • • •■ 
The algorithm that computes the triples St„, Wt„} is the following: 

1. We denote St^ = s and define the set Tg of target states connected to s: Tg — {s', Ts's > 
0}. We also define the total width Tg — J2s'eT, ^s's- 

2. Extract r > with probability density PsiT) — Tge'^''^. In other words, r = — log^ 
with ^ uniformly distributed in [0, 1]. 

3. Extract a new state s' e Tg with probability Ps' — Ts's/^s- 

4. Define t^+i = tn + T, St„^, = s' and Wt„^, = Wt„ ■ e"'^^^. 

In this sector there is no systematic error due to a finite evolution time. The semigroup 
dynamics is in fact reproduced exactly by the above stochastic process (st, Wt). 

About Importance Sampling, wc remark that in the discrete case, the inclusion of a trial 
wave function amounts to the redefinition 

Hs,s = ^^>Hg,g^, (5.14) 

s 

where = (s|\E'q) are the components of the trial ground state \^q). The new Hamil- 
tonian H is not symmetric, but the above formulas works as well with no need for further 
modifications: actually, they have been derived without requiring any symmetry condition 

Some final comments are in order about the choice of the basis for the fermionic states. As 
we mentioned in the general discussion, we want to have zero or negative off-diagonal matrix 
elements of H. The simplest choice amounts to consider eigenstates |n) of the occupation 
numbers xlxi with a relative phase fixed by the natural choice 

I^) = n(x^ri0)' X.|0) = 0. (5.15) 

i=l 
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This does not guarantee that the above sign problems are absent. In fact, in weak-couphng 
perturbation theory, the choice of periodic boundary conditions does not break supersym- 
metry when Lmod4 = as can be checked, e.g., in the free model. However, under this 
condition, there is an even number of fermions, L/2, in the ground state and a sign prob- 
lem arises due to boundary crossings of a fermion, since such a transition involves an odd 
number of fermion exchanges. To avoid such a difficulty, we shall adopt open boundary 
conditions. With this choice, L needs just to be even to assure a supersymmetric weak 
couphng ground state. Also, we shall restrict to the case L mod 4 = 2 and to the sector with 
L/2 fermions that contains a non-degenerate ground state, with zero energy at all orders in 
a weak-coupling expansion. 



D. Algorithm for the full model 

To study the full Hamiltonian of the Wess Zumino model, the simplest attitude is to 
perform an approximate splitting of the bosonic and fermionic sectors. For instance, with 
second order precision, we can write 

exp {-sH) = exp sHb^ exp {-eHp) exp (^-i eHb^ + 0{e^), (5.16) 

and consider separately the evolution related to Hb and Hp freezing the fermionic or bosonic 
fields respectively. In the end, an extrapolation to the £ — > limit must be performed. 
Eq. (5.16) has been approximated to the same order as Eq. (5.9); if necessary, both can be 
improved. 



E. Variance control 



A straightforward implementation of the above formulae fails because of a numerical 
instability: the variance of the walker weights Wt computed over the walkers ensemble 
grows exponentially with t and forbids the projection onto the ground state [27]. A good 
trial wave function can certainly reduce the growth rate, but the problem disappears only 
in the ideal case when the trial wave function is exact. To bypass this problem, some kind 
of branching procedure must be applied in order to delete trajectories (walkers) with low 
weight and replicate those with larger weight. 
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In practice, we introduce an ensemble, i.e., a collection of K independent walkers s^") 
each one carrying its own weight W^'^^: 

E = \<n<K}. (5.17) 

When the variance of the weights in the ensemble becomes too large, £, is transformed in 
a new ensemble 8.' that reproduces the same expectation values (at least in the K ^ oo 
limit) and has a smaller variance. We adopted the branching procedure of Ref. [28]: for each 
walker s^") we compute a multiplicity 

where ^ is a random number uniformly distributed in [0,1], K is the desired number of 
walkers, and \x\ is the maximum integer not greater than x; the new ensemble E' contains 
M^") copies of each configuration s*^") in E and all the weights are set to 1; the actual number 
of walkers K will oscillate around K. This procedure has the advantage that there is no 
harmful effect from its repeated application; therefore we apply it at each Monte Carlo 
iteration, after all the fields have been updated. 

F. Choice and dynamical optimization of the trial wave function 

About the choice of the trial wave function, we propose the factorized form 

1^^) = e^^(^)+^^(^'>^'>^^)|^o)B ® |*o)f, (5.19) 
where |^'o)b ® |^o)f is the ground state of the free model given explicitely in App. (D) and 

= E E = E(-l)" (xiXn - ^ ) E C^k^n- 

n k=l n ^ ^' k=l 

Since the trial wave function is a modification of the free ground-state wave function, we 
expect that importance sampling will improve as we approach the continuum limit. 

The degrees (Ib and dp must be chosen carefully to achieve a balance between the accuracy 
of the trial wave function on one hand, and convergence of the adaptive determination of 
the parameters a and computation time on the other hand. We chose ds — dp — 4, except 
for situations very close to the continuum limit, e.g., V — \2^'^ with A2 < 0.2, for which 
we chose ds — dp — 2 (cf. Sect. VXD). 
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The trial wave function should of course respect the symmetries of the model; a 
symmetry possessed by the model for specific forms of V is very helpful to reduce the number 
of parameters that we must optimize. For odd V ^ the model enjoys the exact symmetry 
—^n, and therefore odd aj^ and can be set to zero. For even V, the model enjoys 
the approximate symmetry — V^m Xn ^ xli (i^ broken by irrelevant terms and by 

boundary terms), and we verified that odd and even can be set to zero. 

Let us denote by a = {a^,a^} the collective set of free parameters appearing in the 
trial wave function. A possible approach consists in performing simulations of moderate 
size at fixed ot in order to optimize their choice. However, as shown in [29], the trial wave 
function |^(^) can also be optimized dynamically within Monte Carlo evolution with a better 
performance of the full procedure. 

The idea is again simple: consider the ground state energy as a typical observable; for a 
given choice of oc, after N Monte Carlo steps, a simulation with an average population of 
K walkers furnishes a biased estimator Eo{N, K,cx.). If we denote by (•) the average with 
respect to Monte Carlo reahzations, Eq is a random variable such that 

lim {Eo{N, K, a)) ^Eo + SEoicx, K), (5.20) 

iV— >oo 

where dEQ{a., K) depends on a, but vanishes as i^T — > oo. Besides, the size of the fluctuations 
is measured by 

Var ^o(A^, K, a) (5.21) 

In the K ^ 00 limit, {Eq) is exact and independent on ol. The constant C2(i^, a) is related 
to the fluctuations of the effective potential V and is strongly dependent on ol. The problem 
of finding the optimal trial wave function can be translated in the minimization of C2{K, ol) 
with respect to ol. 

The algorithm we propose performs this task by interlacing a Stochastic Gradient steepest 
descent with the Monte Carlo evolution of the walkers ensemble. At each Monte Carlo step, 
we update a„ — > a„+i according to the simple law 

otn+i = ttn - ?7nVaVar£:„ V (5.22) 

where £n is the ensemble at step n and {77} is a suitable sequence, asymptotically vanishing; 
to keep things simple, we use the same 77 for all components 0;^ of a, although in principle 
we could use a different 77^ for each a^. 
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A non-linear feedback is thus established between the trial wave function and the evolution 
of the walkers. The convergence of the method can not be easily investigated by analytical 
methods and explicit numerical simulations are required to understand the stability of the 
algorithm. In [29], examples of applications of the method with purely bosonic or fermionic 
degrees of freedom can be found. Here, we apply the method for the first time to a model 
with both kinds of fields. 

In practice, the choice of the initial values of a. is important: it is clear that, if we have a 
good guess of the optimal values (e.g., from runs at the same V but for smaller values of L or 
K), starting from them makes the convergence much faster. We also noticed that, starting 
e.g. from olq = 0, the steepest descent at times fails and a„ oscillates wildly; this never 
happens if most of the starting values have at least the right sign and order of magnitude. 

We found it useful to determine r] dynamically as well. The basic idea is that we wish to 
decrease rj when all the ak have reached the optimal value and they are oscillating at random 
around it: in this situation, a smaller 77 means less noise on a. On the other hand, we wish 
to increase rj when one or more ak is drifting: a larger r) now means a faster approach toward 
the optimal value. To monitor the trend of a, we use the quantity 

T — maxTjfc, Tjfc = o, [p. 16) 

k Vk 

where = ni — no is the number of iterations in the interval of Monte Carlo iterations 
(no,ni] we are considering, Vk is the variance of ak in the interval, and is the slope 
of the linear least-squares fit to a^^n vs. n. is invariant under translations and scale 
transformations; it is positive if ak is drifting and it is negative if ak is oscillating. 

A typical example of the implementation of the a and 77 dynamics is shown in Fig. 4; 77 is 
initiahzed to 10~^; after each interval oi N — 5000 iterations, if r > 0.15, 77 is multiphed by 
V^TO; if r < 0, 77 is divided by \/T0; ifO<r<0.15,77is unchanged; finally, 77 is restricted 
to the interval [10~^,10~^]. The parameters of the r] dynamics given here were obtained 
empirically. 

G. Observable measurement 

We measure the vacuum expectation values of (fnfm, xliXn, and of 

T = i j^iVn+l - Vn-l) V{ipn), Y, = {Q , <pl^2,n} ■ (5.24) 

n=l n=l 
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FIG. 4: a and 7] dynamics from a run at V = 0.5 — 0.55, L = 34 and K = 100. 



Note that, with our choice of boundary conditions, we don't have translation invariance and, 
e.g., will depend on i; however, the dependence is sizable only within a few correlation 
lengths from the border; therefore we typically average site-dependent quantities excluding 
sites closer than a suitable Lmin from the border; in the case of {{pn^Pm), we average over all 
pairs with fixed distance r = \n — m\, excluding the cases when n or m is closer than Lmin 
from the border. 

The ground-state energy is measured simply by averaging the measured values of Eq over 
the ensembles S{t), discarding a suitable thermalization interval {0,to): 



1 



tl-tot=to+i l^i=iWi^t \Si,t\Si,t) 

cf. Eq. (5.2). The vacuum expectation value of a generic observable is computed implement- 
ing the forward-walking formula (5.3) as 



E., 



\Si,t\H\si^t) 



i,t 



(5.25) 



In principle, in Eq. (5.26) we must take the At — > 00 limit, but in practice a moderate value 
like At = 500 is sufficient. A typical examples of At dependence is shown in Fig. 5; it should 



(5.26) 
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FIG. 5: The central charge T (cf. Eq. (5.24)) vs. At from runs atV = 0.5 and L = 34. 



be noticed that the error bars grow with At but very slowly. 



VI. NUMERICAL RESULTS 



A. Review of previous lattice studies 

The class of models that we study in this paper has been previously considered in [18] 
with a Monte Carlo approach that determine the ground state energy by using 

Tr(i7e-^-^) 



En = lim 



(6.1) 



Tr(e-/5^^) ' 

and working numerically with a large fixed jS. This is in the spirit of the usual Lagrangian 
algorithms to be compared with the Green Function Monte Carlo method where j3 can be 
identified with the simulation time and is thus taken to infinity by the very nature of the 
algorithm. 

The analysis of [18] is performed on 12 x 100 lattice, hence with a rather coarse spatial 
mesh. In the model with V{(p) = X^ip^ supersymmetry appears to be unbroken in full 
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FIG. 6: The ground-state energy density Eq/L vs. 1/K V = cp^, L = 22, with statistics of 1 M 
iterations for K < 5000, 500 k iterations at = 5000, and 300 k iterations at = 10000. 

agreement with our analysis. In the quadratic model with V{ip) = Xi^"^ + Aq, the authors 
of Ref. [18] find rather strong signals for supersymmetry breaking with Aq bigger that the 
critical value Aq — —0.5 and have numerical results showing a very small ground state energy 
for Aq < —0.5. No discussion of the continuum hmit is attempted. 

B. Odd V 

As an example of odd prepotential, we study the case V = ^p^ . We plot the ground-state 
energy vs. K in Fig. 6 and the Ward identity Yi vs. K in Fig. 7. Both give a very convincing 
evidence for unbroken SUSY; all the other Ward identities are consistent with zero, but more 
noisy. It should be noticed that the bosonic and fermionic contribution to are ~ ±7.4: 
we are observing a cancellation of four orders of magnitude. 
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FIG. 7: The Ward identity Yi vs. 1/K at V = ip^, L = 22. 



C. Super symmetry breaking 



We now turn to the more interesting case of even prepotential, investigating the case 
V = X2f'^ + Aq. We remind that in this case the model enjoys an approximate Z2 symmetry 
iPn — V^n, Xn ^ Xn- fixed A2, we may expect (in the L — > cxd hmit) a phase transition 
at Ao = Ao'^''(A2), separating a phase of broken SUSY and unbroken Z2 (high Aq) from a phase 
of unbroken SUSY and broken Z2 (low Aq). We investigated in details the case A2 = 0.5. 

The usual technique for the study of a phase transition is the crossing method, applied 
to the Binder cumulant [30] 



B 



1 



2 (M2)2 J ' ^^-^^ 
in our case, the choice of a sensible definition for the magnetization M is nontrivial, since 
our model is neither ferromagnetic nor antiferronagnetic, and it doesn't enjoy translation 
symmetry. We tried out several definitions, before making the choice 



even (odd) 



L - 2Lr 



E 



(6.3) 



^mm gygn (oj^) i=l+Lynin 

where, typically, Lmin = 6; Meven and Modd are perfectly equivalent, and the values of B we 
quote in the following are the average of -Seven and -Bodd- 
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The crossing method consists in plotting B vs. Aq for several values of L; the crossing 
point Aq''(Li, L2), determined by the condition 

B{Xo, Li) = B{Xq, L2), 

is an estimator of Aq'^'' [30]; the convergence is dominated by the critical exponent u of the 
correlation length and by the critical exponent lu of the leading corrections to scaling (cf. 
Ref. [31]): 

A-(Li, L,) = Af) + o(L^-'/^ Lr-'^i; 

we expect the phase transition we are studying to be in the Ising universality class, for 
which v = 1 and uu = 2, and therefore we expect fast convergence of Aq'^ to Aq'^^ The results, 
plotted in Fig. 8, indicate A^ ^ = -0.48 ± 0.01. 

It is possible to study the phase transition by looking at the connected correlation function 
Gd — {(Pnfm)c, averaged over all n, m pairs with \m — n\ — d, excluding pairs for which m 
or n is closer to the border than (typically) 6. In our staggered formulation, even and odd 
d may correspond to different physical channels. 

Gd is fitted to the form exp[— oi — a2d + 03/ {d + 10)], separately for even and odd d; the 
best fits give a good if we remove the smallest distances, typically 0? < 3 for the odd 
channel and d < 4 for the even channel. In the broken phase, we have small but nonzero 
02, and we observe equivalence of the even- and odd-d channels; an example is shown in 
Fig. 9. In the unbroken phase, 02 is larger, and the even- and odd-d channels are somewhat 
different; an example is shown in Fig. 10. The difference between the two phases is apparent, 
e.g., in the plot of 02 vs. Aq, presented in Fig. 11; the data presented here would indicate 
Af^ = -0.48 ±0.01. 

An alternative window to the phase transition is offered by the optimized values of the 
parameters of bosonic part of the trial wave function, which should be related to the effective 
potential Ves of the bosonic field: 

we verified that ctf < 0, as required by stability; a negative value for therefore indicates 
unbroken Z2 symmetry, while a positive ctf indicates spontaneous breaking of Z2. 

The numerical values of are shown in Fig. 12. It is clear that, especially in the broken 
phase, statistical errors are underestimated. The abovementioned scenario is qualitatively 
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FIG. 8: The Binder cumulant B vs. Aq- 



(c) 

confirmed, but the data yields Ag ~ —0.40, which is very far from the more traditional 
estimates obtained above. 

Finally, to investigate the supersymmetry properties of each phase, we analyze S, the 
ground-state energy density extrapolated to infinite K and L. We fit Eq/L to the form 



Eo _ ^ , a + bL 



(6.4) 
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FIG. 9: The connected correlation function Gd for V = 0.5 99^ — 0.5. 




FIG. 10: The connected correlation function Gd for V = 0.5 99^ — 0.38. 
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FIG. 11: The effective mass 02 of Gd vs. Aq for L = 58 and K = 500; for Aq < —0.51 the error on 
02 is very large. 
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FIG. 12: The optimization parameter a^; it is insensitive to L and K. 
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FIG. 13: The ground-state energy density Eq/L vs. Aq- The sohd hne is a fit to the form Eq/L = 



X^/#d.o.f from 1 to 2; the errors on the fit parameters are defined by the values giving an 



increase of by 1; if > #d.o.f we multiply them by the "scale factor" S = yx^/^d.o.f. 



(c) 

We present a plot of S vs. Aq in Fig. 13; these data give Xq ~ —0.53, with a rather large 
uncertainty. 

D. Continuum limit 

We wish to investigate if the pattern established in Sec. VI C extends to the continuum 
limit. 

We study the trajectory 



corresponding to a 1-loop RG trajectory, cf. Eq. (4.6); the effect of Aq is small in the range 
we considered, therefore we expect this to be a reasonable approximation to a true RG 
trajectory. 

We estimate the correlation length from the exponential decay of the connected correla- 
tion function Ga = {fn'^m)c averaged over all n,m pairs with \m — n\ = d, excluding pairs 





Ao = ln(4A2), 

ZTT 



(6.5) 
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for which m or n is closer to the border than (typically) 8. In our staggered formulation, 
even and odd d correspond to different physical channels. 

We performed runs for values of A2 spaced by a factor of -s/2, with a statistics of 4x10^ 
iterations. In Figs. 14 and 15 wc show the plots of the (/? correlation for the case V = 
0.35 (f^ + 0.02. It is very difficult to extract a correlation length from the even-ci channel, 
presumably because (p has a very small overlap with the lightest state of the channel, and 
the value 1/^ — 0.20 ± 0.03 quoted in Fig. 14 should be considered tentative. The odd-d 
channel is much cleaner, and it is possible to estimate ^ with a good precision. 

For the other values of A2, the situation is similar but with shghtly larger errors. The 
measured values of ^odd follow very weU the naive scaling behavior 



The entire range 0.088 < A2 < 0.35 seems to be in the scaling region, with A2 = 0.5 a 
borderline shown in Fig. 16. The values of ^even have very large errors, and it is 

hard to draw any conclusion from them. 

The Green Function Monte Carlo algorithm gives a very accurate measurement of the 
ground-state energy Eq; to give a feeling of the precision reached, wc quote the results for 
the smallest and the largest value of A2 we considered, along the trajectory (6.5): 



Eq show a sizable dependence on K and L, and it is therefore necessary to perform an 
extrapolation to L ^ 00 and K ^ 00; we fitted the energy density to the form 



which gives a good x^- ^ remains constant within errors, with a value u ~ 0.75, i.e., the 
algorithm performs well as we approach the continuum limit. 

The "scaling" plot of the energy density Eq/L is shown in Fig. 17. It seems to behave as 
A2'^, while naive scaling would predict Eq/L oc A|. 

The nonzero value of Eq/L (disregarding this puzzling exponent) and the lack of any 
signal for a breakdown of parity show that the trajectory we are considering belongs to the 
phase with broken supersymmetry and unbroken Z2 symmetry. 



e oc I/A2. 



Eo(A2=0.044,L=46,ir=200) = (1.28 ± 0.01) x 10 



-3. 



£;o(A2=0.5,L=46,is:=200) = (69.44 ± 0.05) x 10 



-3 
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d 

FIG. 14: The connected correlation function at even distance for V = 0.353553 (/j^ + 0.019502; 
the curve and value of 1/^ quoted are the result of an exponential fit for 10 < d < 18 to the L = 46, 
K = 200 data. 

VII. CONCLUSION 

In this paper, we investigated a class of two-dimensional = 1 Wess-Zumino models 
by non-perturbative lattice Hamiltonian techniques. The key property of the formulation 
is the exact preservation of a SUSY subalgebra at finite lattice spacing. Our main tool are 
numerical simulations using the Green Function Monte Carlo algorithm; we also performed 
strong-coupling expansions. 

All our results for the model with cubic prepotential indicate unbroken supersymmetry. 
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the curve and value of 1/^ quoted are the result of an exponential fit for 3 < d < 15 to the L = 46, 
K = 200 data. 

We studied dynamical supersymmetry breaking in the model with quadratic prepotential 
V = A2V?^ + Aq, performing numerical simulations along a line of constant A2. We confirm 
the existence of two phases: a phase of broken SUSY and unbroken Z2 at high Aq and a 
phase of unbroken SUSY and broken Z2 at low Aq, separated by a single phase transition. 

We studied the approach to the continuum limit in the model with quadratic prepoten- 
tial performing numerical simulations along a 1-loop RG trajectory in the phase of broken 
supersymmetry; we measured the correlation length of the bosonic field (in the odd-distance 
channel), which is found to scale with the expected exponent; on the other hand, the ground- 
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FIG. 16: The correlation length at odd distance ^odd- The dashed curve is the result of a scaling 
fit (with fixed exponent). 
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FIG. 17: The ground-state energy density Eq/L, extrapolated to L ^ cxd and K oo. The dashed 
curve is the result of a two-parameter fit, while the solid curve shows the naive scaling behavior. 
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state energy density scales with an exponent clearly different from the expected exponent. 

In many instances, the simulation algorithm suffers from slow convergence in the number 
of walkers K. 
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APPENDIX A: CHECK OF UNBROKEN SUSY FOR V{ip) = X-s_ip + Aq 



If the potential V{(p) is a linear function of the field if, then the ground state can be 
found explicitely. With a field translation we can set Aq = 0. The model Hamiltonian is 
Hb + Hp where we recall that 



n=l 
L 



2^" + 2 I 2 + 



1 



n=l L ^ 

Thus, in the bosonic sector, the Hamiltonian can be written 

1 1 

Hb^ -YpI + I^Y fnVnmVm: 
n nm 

with 



(Al) 
(A2) 



(A3) 



B 



(A4) 



\l + 1/4 n = m and n = 1,L 
Xl + 1/2 n — m and 1 <n < L 
-1/4 |n - m| = 2 

In the fermionic sector, the Hamiltonian can be written in terms of canonical Fermi annihi 
lation and creation operators as 



(A5) 



with 



(-1)"A2 n = m 



Xl + 1/2 n — m and 1 <n < L 
-1/2 |n-m| = l 



(A6) 
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If we denote by {{co^y} and {co^} the sorted eigenvalues of and V^, then we find that 
the ground state has actually zero energy 

L L/2 
^ n=l n=l 

This can be proved in the spirit of SUSY without computing explicitely the eigenvalues. In 
fact, we can check that {V^Y is the matrix apart from a wrong sign in the diagonals 
|n — m| = 2. This can be repaired by changing sign ip ^ —ip on the sites with ?T,mod4 = 1, 2. 
Taking into account the particle-hole symmetry of Hp, we have thus proved that the spectra 
of criV^) and have the general form 



aiy^) = {-Xi,Xi, -X2, X2, . . .} 
a{V'')^{xl,xl,xl,xl,...} 



(A8) 
(A9) 



with full cancellation between the lowest L/2 fermionic values and one half of the square 
root of the bosonic ones as in Eq. ( A7) . 

APPENDIX B: STRONG-COUPLING EXPANSION OF {ipk) AND {ipkipi)c 
1. {<^k) 



Let us define 



^0 



The vacuum expectation value of the field </? is 



(1)| 



^{-lf{l-2{m'i\nk\^\,")). 



(Bl) 



(B2) 



The expectation value of the occupation number can be computed by going to the basis {a} 
and is 



L/2 



A straightforward calculation gives 



(B3) 



1 

2L 



cos 



1 + L + 



TT 

2L 



{2k{L + l) - 1 



TT 



(B4) 
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and therefore 



cos 



TT 



— (2/c(L + l)-l 



TT 



Sin— (2A;-1) 



> . 



(B5) 



It is interesting to consider the hmit L — > oo of this expression after a rescahng A; — > xL 
where < x < 1. The result is 



^ ^±1 + cot(7ra;)) + ^ 



1 



2^/2eo{L^ ' ' 2L2sin2(7rx) 

where the sign is +1 for even k = xL and —1 for odd k. 



O 



L3 



(B6) 



2. {(Pk<Pl)c 



Let us denote briefly 



and 



The 2-point correlation is, lor k ^l, 

Going to the {a} basis, we immediately obtain 



{nkni)c= '"kvf- I] 



k ^/ > 



l<A<Lj2 



L/2+l<B<L 



and, for k ^ I, 



(B7) 

(B8) 

(B9) 



(BIO) 



(Bll) 



The two sums over eigenvalues can be evaluated analytically thanks to the simple form of 

(v) 

the eigenvectors v} . The explicit result is 



T L/2 J . 

p=l+L/2 



(B12) 
(B13) 
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where 



n even, m even : (—1) 2 



n even, m odd : (— 1)"^2' 



TT 



1 + cot -^{tn + n — 1) 



TT 



1 + cot — (n — mj 



n odd, m even : (—1) 



TT 



1 + cot 777 — 



(B14) 



n odd, m odd : (— 1)"2'' 



TT 



-1 + cot — (m + n - 1) 



that can be used to compute the connected correlation on a finite lattice. 

It is interesting to note that the limit L ^ 00 can be taken without rescaling n and m. 
For instance, we have 

1 

4(^)2 



lim (^'fi'Vl V^nl^O 



n even : 



n odd 



(n - 1)2 
1 



(B15) 



APPENDIX C: SECOND-ORDER EXPANSION FOR Eq, EVEN q 

The general formula for the second order contribution to the ground state energy is 

E2 = £"2,1 + -£'2, 2 
-E'2,2 



(^W|i/(2)|^')(^'|i7(2)|^i: 



E 



Eq-E' 



The states |\E'') are excited states of the form 



l*') = CW---CWK'---^^) 

where ki-\ ki — v > ^ (integer) and ui — (—1)"'+'. For such a state we have 



(CI) 
(C2) 

(C3) 



A first important remark is that if^^) (.9,^ restricted to its H^^^ part. In fact, the fermionic 
part of i?(2) can be written as Hf^^ plus terms that are only functions of the occupation 
numbers. These operators acting on give states that are orthogonal to the subspace 

of excited unperturbed states |^''). 
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The first term in E2 is simple and can be computed straightforwardly exploiting the fact 
that the expectation value of cpf over |^(^)) does not depend on I: 

^2,1 = - \ E V^m2|*f ^> = - Im'coiL) (C4) 



where 

/oo 
-00 

and 

Co(L) = E(l - (^/> - (^'+2) + Mriini+2)) (C5) 



About -^2,2, it can be computed by considering states with an excited single-site wave func- 
tion in one or two distinct sites. Summing the two contributions we find 

s>0,t>0 ^ l.^ L 

+ C3(L)$«$fVWy/^)| (C6) 
The symbols appearing in the above equations are defined as follows 

= (V^o'=l^(<^)lV'o'=) = V^Torio (C7) 
Vj'^ = -^iV^o + i-iyV^ii^M), (C8) 

The functions co,i,2,3(-^^) can be fitted with a few powers of L and the numerical result is 

co(L) =L- 1.495 - + • • • (CIO) 

= -0.540- If? + H?! + ... (CU) 
c.(L) ^ iL- 0,540 + !^ + H^^... (C12) 

C3(L, = -|l-^ + ... (C13, 
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The full expression for E2 has thus the following simpler asymptotic expression 



, 2 



+ 



- \ . (C14) 



s>0,t>0 

As an example, we find for the model with V{ip) — ip^: 

^,,2 , 0.160 0.0481 , ^,,2 N ^o... 0.160 0.0488 

£;(A^ 22) = 0.2811 - + and ^(A^, 00) = 0.2811 - + 

If we are interested in a comparison with an actual simulation, some trivial rescaling is 

necessary. For instance, in the case of a purely quadratic V{ip) = A2 we must compare 
the expansion and the results for X^^^E^jL and identify A = A2^^ as the strong-coupling 
expansion parameter. 

APPENDIX D: FACTORIZED WAVE FUNCTION FOR THE FREE MODEL 

In the free case F = 0, with even L, we have 

//(•^^ = //b + Ep, (Dl) 
where ((/?o = V'l+i = 0, the same for fermions) 

= E [IpI + li^n+l - Vn-l)'} , (D2) 
1 ^ 

Hf^ -7^J2 (XnXn+l + xl+lXn) ■ (D3) 

^ n=l 

The two Hamiltonians Hb and Hp are decoupled and the ground state takes the factorized 
form 

|*(o)) = |*g))®|*g)). (D4) 

We now determine 1^^''^;') separately, assuming L mod 4 = 2 in order to have a unique ground 
state in the decoupled = model. 
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1. Bosonic sector 



Let us write Hb in the form 



1^ o 1 ^ 



^J2Pn+2 51 Vnm'Pn'Pm- (D5) 

n=l n,m=l 

Let A| be the eigenvalues of the L x L matrix V and let z^^ be the corresponding (real) 
eigenvectors satisfying 

z^'^ ■ z^'^ = 5,,^. (D6) 

The ground state is 

=exp [-^ 5^ i^nmV'nV'm) , (D7) 
Y n,m=l j 

where 

Rum = E A,4'^4'^ (D8) 
fe=l 

The explicit form of A^; is 

Afe = sin(|^0, ^=Li(A; + l)J (^ = l,...,iL); (D9) 
the dimension of each eigenspace is 2. The ground state energy of Hb is 



^o,B = E|A. = E«iM-^^j. (DIO) 



We adopt for the two orthogonal eigenvectors the choice 



2 (n even) 

V-f' + l 1^ sm (^ ^_^ ^ -(L - n + 1)J (n odd) 



(2i) ^ ^ 



. /2^-l7r \ 

sm — n \n even) 

VL + l 2 ; ^ ^ (D12) 

(n odd) 



2. Fermionic sector 

In the free case, the Hamiltonian in the fermionic sector is diagonalized by the orthogonal 

change of basis 

' n=l 
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and takes the form 



Hp ^-2^ cos Y<an, 



n=l 



Hence the one-fermion energies are 



■ cos ■ 



n— 1, . . . , L. 



The fermionic component of the (supersymmetric) ground state is 

L/2 



^f) = n 4|0); 



n=l 



the ground state energy of Hp is 



L/2 



£;o,^ = -gcos^; 



(D14) 



(D15) 



(D16) 



(D17) 



we can easily check that 



Eq^f — —Eq^b — 



TT 



sm 



(D18) 



V 2(L + 1)/ 



and therefore Eq — Eq b + Eq p — 0- 
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